package xj.graph2d.util;

import java.awt.image.BufferedImage;
import java.awt.image.ConvolveOp;
import java.awt.image.Kernel;

public class GaussianBlur {

  public GaussianBlur(double stdDev) {
    this(stdDev, stdDev);
  }

  public GaussianBlur(double stdDevX, double stdDevY) {
    this.stdDevX = stdDevX;
    this.stdDevY = stdDevY;
    xinset = surroundPixels(stdDevX, highQuality);
    yinset = surroundPixels(stdDevY, highQuality);

    if ((xinset != 0) && ((stdDevX < 2) || highQuality)) {
      convOp[0] = new ConvolveOp(makeQualityKernelX(xinset * 2 + 1));
    } else {
      dX = (int) Math.floor(DSQRT2PI * stdDevX + 0.5f);
    }

    if ((yinset != 0) && ((stdDevY < 2) || highQuality)) {
      convOp[1] = new ConvolveOp(makeQualityKernelY(yinset * 2 + 1));
    } else {
      dY = (int) Math.floor(DSQRT2PI * stdDevY + 0.5f);
    }
  }

  protected double stdDevX, stdDevY;

  protected ConvolveOp[] convOp = new ConvolveOp[2];

  boolean highQuality = true;

  int xinset, yinset;

  int dX, dY;

  /**
         * Constant: sqrt(2*PI)
         */
  static final float SQRT2PI = (float) Math.sqrt(2 * Math.PI);

  /**
         * Constant: 3*sqrt(2*PI)/4
         */
  static final float DSQRT2PI = SQRT2PI * 3f / 4f;

  /**
         * Constant: precision used in computation of the Kernel radius
         */
  static final float precision = 0.499f;

  public static int surroundPixels(double stdDev, boolean highQuality) {
    if ((stdDev < 2 || highQuality)) {
      // Start with 1/2 the zero box enery.
      float areaSum = (float) (0.5 / (stdDev * SQRT2PI));
      int i = 0;
      while (areaSum < precision) {
	areaSum += (float) (Math.pow(Math.E, -i * i / (2 * stdDev * stdDev)) / (stdDev * SQRT2PI));
	i++;
      }

      return i;
    }

    // compute d
    int diam = (int) Math.floor(DSQRT2PI * stdDev + 0.5f);
    if (diam % 2 == 0) {
      return diam - 1 + diam / 2; // even case
    } else {
      return diam - 2 + diam / 2; // Odd case
    }
  }

  private float[] computeQualityKernelData(int len, double stdDev) {
    final float kernelData[] = new float[len];

    int mid = len / 2;
    float sum = 0; // Used to normalise the kernel
    for (int i = 0; i < len; i++) {
      kernelData[i] = (float) (Math.pow(Math.E, -(i - mid) * (i - mid)
	  / (2 * stdDev * stdDev)) / (SQRT2PI * stdDev));
      sum += kernelData[i];
    }

    // Normalise: make elements sum to 1
    for (int i = 0; i < len; i++) {
      kernelData[i] /= sum;
    }

    return kernelData;
  }

  private Kernel makeQualityKernelX(int len) {
    return new Kernel(len, 1, computeQualityKernelData(len, stdDevX));
  }

  private Kernel makeQualityKernelY(int len) {
    return new Kernel(1, len, computeQualityKernelData(len, stdDevY));
  }

  public BufferedImage filter(BufferedImage src) {
    BufferedImage des = src;
    if (src != null) {
      if (convOp[0] != null) {
	des = convOp[0].filter(src, null);
      }
      if (convOp[1] != null) {
	des = convOp[1].filter(des, null);
      }
    }
    return des;
  }

}
